Duality, thermodynamics, and the linear programming problem 
in constraint-based models of metabolism 
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It is shown that the dual to the linear programming problem that arises in constraint-based models 
of metabolism can be given a thermodynamic interpretation in which the shadow prices are chemical 
potential analogues, and the objective is to minimise free energy consumption given a free energy 
drain corresponding to growth. The interpretation is distinct from conventional non-equilibrium 
thermodynamics, although it does satisfy a minimum entropy production principle. It can be used 
to motivate extensions of constraint-based modelling, for example to microbial ecosystems. 
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In biology, the metabolism of an organism provides 
energy and raw materials for maintenance and growth. 
As such, an interesting and important question concerns 
the application of thermodynamics to metabolic reaction 
networks [H, 0, H, Hj] . For example, Prigogine and Wiame 
suggested a long time ago that an organism's metabolism 
might be governed by a minimum entropy production 
(MEP) principle [|. From the physical point of view, 
a metabolic reaction network is an excellent example of 
a system in a non-equilibrium steady state, since one 
can usually assume that the metabolite concentrations 
are unchanging after a short transient relaxation period. 
The appropriate generalisation of thermodynamics and 
statistical mechanics to non-equilibrium steady-states is 
a large field which continues to attract attention to 
this present day 0]. In this Letter, we show that a novel 
thermodynamic interpretation can be given to the dual 
linear programming problem which arises in constraint- 
based models of metabolism. The resulting interpreta- 
tion is rigorously defined, and uniquely determined by 
the mathematics. It is closely analogous to, but distinctly 
different from, conventional non-equilibrium thermody- 
namics. We also show that it satisfies an MEP principle 
similar to that proposed by Prigogine and Wiame. 

Constraint-based modelling (CBM) of metabolic net- 
works has been pioneered by Palsson and co-workers Q. 
In a typical application, described in more detail below, 
the steady-state assumption is combined with a target 
function to make a linear optimisation or linear program- 
ming (LP) problem. The LP variables are the fluxes 
through the various reactions that comprise the network, 
and the LP constraints arise from basic considerations of 
stoichiometry and from the reversibility or otherwise of 
the reactions. The LP objective function is biologically 
motivated, for example a 'growth' reaction is commonly 
inserted, and the target is to maximise flux through this 
reaction to correspond to maximal growth rate. CBM has 
been applied to microorganisms from all three domains 
of life [3, III 11 1 , and has been remarkably successful in 
predicting phenotypic behaviour T3, IH, 3- 
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Mathematically, every LP problem has a unique dual 



FIG. 1: Schematic metabolic network for a prokaryote like E. 
coli, showing intracellular metabolites (open circles), extra- 
cellular metabolites (hatched circles), internal and exchange 
reactions (arrows), and a growth reaction (dashed arrow). 



[H| • It was in determining the dual to the CBM LP prob- 
lem that we noticed a striking analogy to non-equilibrium 
thermodynamics. Let us start therefore with a general 
discussion of LP duality, before specialising to the case of 
CBM. We recall that the basic or primal LP problem, in 
standard form, is to maximise an objective function z = 
ELi a a x a given J22=i A iaXa = h (i = 1 . . . m, m < n), 
where the x a > are variables, the a a are coefficients, 
Ai a is a matrix, and the bi are constants (we use Greek 
and Roman indices to emphasise that different compo- 
nents live in different vector spaces). The dual problem 
is then to minimise an objective function w ~ X/j=i 
subject to X/i=i ttiAia > tt a , with no restriction on the 
sign of dual variables Tti. The LP strong duality theorem 
guarantees that max z = min u>, provided both problems 
have optimal solutions. In addition, at optimality, 'com- 
plementary slackness' (CS) conditions hold. To formulate 
these, first define the 'slack' in the inequalities in the dual 
problem to be y a — YliLi n i-^-ia ~ a a- The CS conditions 
state that the inequalities x a y a > are saturated (i. e. 
— 0) at optimality, and only at optimality. 

In many applications of LP, the dual problem can be 
given an economic interpretation, which has led to the 
dual variables being generically known as 'shadow prices'. 
We note that shadow prices can be obtained directly from 
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Primal 



Dual 



CS inequalities 



variables 

flux balance / thermodynamics 
reversible reactions 
irreversible reactions 
growth reaction 
limited exchange reactions 
objective function 



fluxes, v a 

y ] a Si a v a = o 

v a unlimited 

v a > 
v gr > 

Vc > ~V™ n 



chemical potentials, \ii 

B a — ^ \ fliSia 

B a =0 
B a < 

B a < -B* < 
B a < 

W = J2a V a iU \B a \/B* 



v a B a < 

Vgr(B a + B*) < 

(v a + v™ in )B a < 



TABLE I: The primal and dual linear programming problems in constraint-based models of metabolism. In the dual objective 
function, the sum is over the limited exchange reactions only. The complementary slackness (CS) inequalities are saturated 
(i. e. — 0) at optimality, where also max z — min w. The shadow prices for the primal problem are given by — \iijB* . 



the solution to the primal problem [la] , so the dual prob- 
lem need never be explicitly formulated. This may be the 
reason why the remarkably simple structure of the dual 
problem in CBM has not been described before. The use 
of shadow prices in CBM was pioneered by Varma and 
Palsson to assess efficiencies in a model of the central 
metabolism of E. coli fltlj ]. 

Now let us turn to the LP problem in CBM. We start 
with the set of chemical rate equations that describe the 
metabolic reaction network, dci/dt = ^2 a Si a v a , where 
the Ci are metabolite concentrations, the v a are reaction 
velocities or fluxes, and Si a is a stoichiometry matrix giv- 
ing the number of moles of the ith metabolite involved in 
the cvth reaction. Making the steady-state assumption, 
the chemical rate equations reduce to a set of flux-balance 
conditions ^ Q Si a v a = 0. At this point, in CBM, at- 
tention shifts from the metabolite concentrations to the 
reaction fluxes. From this point of view, the flux-balance 
conditions become a set of linear constraints on the v a . 
In addition, one usually imposes the 'thermodynamic' 
constraint that v a > if a reaction is irreversible. 

In modern approaches 0] , the reactions in the network 
are elementally and charge-balanced. To make the net- 
work 'do' something, two kinds of imbalanced reactions 
are typically added. The first, as mentioned already, is 
a growth reaction. This reaction drains the endpoints 
of metabolism in the appropriate ratios and represents 
the combined effect of the biochemistry subsequent to 
metabolism. The flux through the growth reaction (the 
growth rate) will be labelled v gI . The second type of im- 
balanced reaction is an 'exchange' reaction, which repre- 
sents the exchange of an extracellular metabolite with the 
environment (the model additionally includes transporter 
reactions which allow extracellular metabolites to enter 
and leave the intracellular environment). The exchange 
reactions enable the uptake of food substrates, trace min- 
erals, dissolved gases, and vitamins; and the discharge 
of metabolic waste products. By convention, a positive 
(negative) flux through an exchange reaction represents 



the discharge (uptake) of the corresponding metabolite. 
Exchange reactions may be reversible, or irreversible if 
discharge only is possible. A special case arises when one 
wishes to represent limited availability, for example of a 
food substrate. In this case the exchange flux is allowed 
to become negative to a limited extent, thus v a > — u™ m 
with i>™ ln > representing a 'cap' on the (negative) re- 
action flux. The value of v™ ln is typically empirically 
determined to agree with experimentally measured up- 
take rates. Fig. [1] shows schematically how the exchange 
reactions and the growth reaction are connected into the 
rest of the metabolic network. For high accuracy work, 
an ATP — > ADP maintenance reaction with a specified 
flux is sometimes included in the model [111 ]. We omit 
this here although it can easily be accommodated with a 
small extension to the formalism. 

The LP problem in CBM is then to find values for the 
fluxes v a which maximise v gI subject to the above con- 
straints. This is summarised in Table [H Usually it is the 
limited availability of substrates through the exchange 
reactions that prevents the problem being unbounded. 
Technically the LP problem is not quite in the standard 
form but it does not take much to make it so. 

We formulate the dual to the above (primal) LP prob- 
lem following the textbook approach described above. 
After some straightforward simplifications, the follow- 
ing picture emerges. Each metabolite has an asso- 
ciated shadow price iii which is unrestricted in sign. 
Each reaction has an associated constraint, of the form 
Y]- TTiSia = for reversible and unlimited exchange reac- 
tions, J2 i TTiSi a > for irreversible and limited exchange 
reactions, and 7TiSi a > 1 for the growth reaction. The 
objective function is w = J2ia ^iSiaV™ 111 where the sum 
is over the limited exchange reactions only. The LP prob- 
lem is to find values for the shadow prices -Ki which min- 
imise this objective function subject to the constraints. 
Note that, although the v™ m appear in the dual objec- 
tive function w, these are numerical constants common 
to both the primal and dual problems. The actual fluxes 
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FIG. 2: Chemical potential analogues (shadow prices) plotted 
as a function of metabolite molecular weight, for the intracel- 
lular metabolites with ^ > in a genome-scale constraint- 
based model of E. colt [IJ. A normalising factor is included 
to make the fa dimensionless. 

v a do not feature in the dual problem. 

We now show that the dual problem admits a strik- 
ing thermodynamic interpretation. The motivation is the 
standard expression for the free energy change in a chem- 
ical reaction, or reaction affinity, A a — ^2 i fJ,iSi a , where 
the Hi are chemical potentials [6| • The similarity between 
this, and the rules for formulating the dual LP problem 
above, make it natural to interpret the shadow prices as 
chemical potential analogues. To aid the interpretation, 
we rescale the dual problem by a factor ~B* < 0, set 
Hi = —iTiB*, and write B a = Y^i^i^ia as the analogue 
of reaction affinity. The resulting thermodynamic formu- 
lation of this rescaled dual LP problem is summarised in 
Table HI We have introduced B a to distinguish our inter- 
pretation from conventional non-equilibrium thermody- 
namics; in general B a ^ A a , as explained in more detail 
below. 

Let us discuss the thermodynamic interpretation in a 
bit more depth. We see that the constraints assert that 
B a = for a reversible reaction, and B a < for an ir- 
reversible reaction. These are precisely in accord with 
equilibrium chemical thermodynamics. In addition we 
interpret the fact that B a < —B* < for the growth re- 
action to mean that a minimum free energy drain equal 
to B* is required for growth. The magnitude of B* sets 
the overall energy scale, and can be arbitrarily chosen. 
Finally, in the rescaled dual LP problem the objective 
is to minimise wB* = ^ Q v™ ln \B a \, in other words a 
weighted sum of the free energy consumption associated 
with the limited exchange reactions. At optimality, one 
has max z — min w, hence the growth rate is easily calcu- 
lated from the solution to the thermodynamic LP prob- 
lem using v gI = £ Q v™*\B a \/B*. 

We now show that the formalism satisfies an MEP 
principle. To derive this, we consider the internal en- 



tropy production due to the chemical transformations, 
TS = — Y^ia v aBa- The sum excludes the exchange re- 
actions since the flux balance condition implies the to- 
tal entropy production ^ Q v a B a = when the sum is 
over all reactions. It is straightforward to show that 
z < TS/B* < w. Thus, at optimality, the entropy 
production is 'pinched' between the two objective func- 
tions. Alternatively, for a fixed growth rate u gr , one has 
v gI B* < TS. Since the minimum value is attained at the 
combined solution of the primal and dual problems, this 
gives the desired MEP principle. 

What does the dual solution look like for a constraint- 
based model of metabolism? To give an example, we 
computed the chemical potential analogues for a genome- 
scale model of the metabolism of E. coli growing aero- 
bically on a glucose 'minimal medium', for which uptake 
of extracellular glucose is the limiting exchange reaction 
[I?! ■ Lack of space precludes a detailed discussion of the 
results, but we find that the vast majority of chemical po- 
tentials are positive and there is a broad distribution over 
several decades of magnitude. An interesting observation 
is that the chemical potentials increase with increasing 
molecular complexity. This is shown in Fig. [2| using 
molecular weight (discounting metal ions) as a stand-in 
for molecular complexity. This correlation arises because 
the chemical potential of a complex molecule is given ap- 
proximately by the sum of the chemical potentials of its 
constituent parts. This in turn follows from the CS con- 
ditions which imply B a — J^. [nSi a — for all reactions 
with a flux v a (see further discussion below). 

Let us discuss our findings in a wider context. Our 
results show that the dual to the CBM LP problem has 
a thermodynamic interpretation in which the dual vari- 
ables are analogous to chemical potentials. The rules 
to formulate this thermodynamic LP problem are sum- 
marised in Table [T] In principle we could strip away the 
CBM 'scaffolding' and let the thermodynamic LP prob- 
lem stand on its own, since LP duality guarantees this 
is equivalent to solving the original (primal) LP prob- 
lem. Such a viewpoint motivates a number of interesting 
questions. 

Firstly, a technical point arises since the primal LP 
problem is often degenerate, in the sense that alternative 
optimal flux distributions exist . This reflects the fact 
that multiple pathways may exist in the metabolism. But 
this is not a serious problem, for the dual problem will be 
similarly degenerate but the strong duality theorem and 
the CS conditions still hold, allowing one to move from a 
solution of the dual problem to a solution of the primal 
problem, and vice versa. 

A more serious discussion point concerns the relation- 
ship to conventional non-equilibrium thermodynamics. 
For the reversible reactions, B a = is a constraint. For 
the irreversible reactions, the CS conditions (Table [TJ 
show that B a — if there is a flux (v a > 0) through 
a reaction, and B a < only if there is no flux (v a = 0) 
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through a reaction (except for the growth reaction where 
we expect v gY > and hence B a — —£?*). This presents 
a sharp contrast to conventional non-equilibrium ther- 
modynamics where a flux through a reaction (v a > 0) is 
associated with a (negative) affinity driving force A a < 0. 
This clearly demonstrates that B a ^ A a , and the ther- 
modynamics described in Table [J is not simply the same 
as conventional non-equilibrium thermodynamics. We 
must therefore regard the rules described in Table Q] as de- 
scribing a novel but tightly constrained thermodynamics 
for the CBM class of problems, derived from the (unique) 
dual to the primal LP problem. Whether the close anal- 
ogy to equilibrium chemical thermodynamics (and the 
unexpected appearance of an MEP principle) is indica- 
tive of deeper principles or not remains a problem for 
future investigation. It would, for example, be an in- 
teresting exercise to compare the E. coli shadow prices 
with what is known about the thermodynamic metabolic 
state of this organism Q • We should emphasise that our 
MEP principle is couched in terms of B a and not A a , 
and therefore our results do not constitute a proof of the 
original proposition of Prigogine and Wiame [B| . 

Another interesting remark is that the choice in the pri- 
mal LP problem to maximise the flux through a growth 
reaction seems to be 'pure biology'. Experiments demon- 
strate that this works well under controlled conditions 
131 ] , and it can be supported by examining population 



helpful discussions and generously providing the E. coli 
model. More details of our derivations, and more exam- 
ples, will be presented in a longer publication currently 
in preparation. 
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dynamics for continuous culture growth in a chemostat 
fl3 ] . Other choices could and perhaps should be made 



in different circumstances [2CJ, [2l|, [22|. In terms of the 
thermodynamic LP problem, this biologically-motivated 
component is translated into the existence of a growth 
reaction with a minimal free energy drain B*. We could 
turn this observation to our advantage, to suggest ex- 
tensions to the CBM approach which are perhaps unob- 
vious in the primal LP problem. Consider for example 
metabolism in a microbial ecosystem, comprising multi- 
ple species which share a pool of common extracellular 
metabolites. The obvious generalisation of the thermody- 
namic LP problem is to include a growth reaction with a 
minimal free energy drain B* for each organism, and seek 
to minimise the free energy consumption of the ecosys- 
tem through the exchange reactions of the extracellular 
metabolites. Of course LP duality means there is a corre- 
sponding primal model (in this case the primal objective 
function becomes a weighted sum of growth rates [Hj])- 
Further exploration of this we leave to future work. 
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